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We study the Josephson current through a resonant level coupled to a vibration mode (local 
Holstein model) in the adiabatic limit of low oscillator frequency. A semiclassical theory is then 
appropriate and allows us to consider the oscillator dynamics within the Born-Oppenheimer approx- 
imation for arbitrary electron-vibration couplings. The resulting Fokker-Planck equation has been 
solved in the most relevant underdamped limit and yields the oscillator distribution function and the 
Josephson current. Remarkably, a transition from single- well to double- well behavior of the effective 
oscillator potential surface is possible and can be tuned by variation of the superconducting phase 
difference. The Josephson current is shown to be only weakly affected by the electron-vibration 
coupling due to strong phonon localization near the bottom of the potential surface. 
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I. INTRODUCTION 

The field of molecular electronics continues to pose 
interesting scientific questions that are also of applied 
relevance. Many aspects of charge transport through 
junctions containing a single molecule have already been 
clarified^ and relatively simple theoretical models 3 - - — 
can often capture the essential physics in such devices; 
see also Ref. [ID for a recent review. To quote just a 
few important experimental works, single-molecule trans- 
port has been studied using different organic molecules^ 
fullerenesrSr— carbon nanotubes ; 23 ' 24 and single hydro- 
gen molecules between Pt leads.— When two supercon- 
ducting (instead of normal-state) electrodes with a phase 
difference (j> are attached to the molecule, the Joseph- 
son effect^ implies that an equilibrium current I(4>) can 
flow through the molecular junction. The impressive ex- 
perimental control over supercurrents through molecular 
junctions achieved recently (see, for instance, Ref. and 
references therein) has been accompanied by first theo- 
retical studies investigating the effects of vibrational or 
conformational molecular modes on the supercurrent. In 
particular, the effect of just one harmonic vibration mode 
coupled to a single-level quantum dot ("local Holstein 
model" ) has been considered in the superconducting ver- 
sion. Analytical results have been obtained via pertur- 
bation theory in the molecule-to-lead tunnel couplings 2 - 
or in the electron-vibration coupling] 29 ' 30 Other works 
have modelled the conformational mode as a two-level 
system i^i 

In this work, we consider the superconducting local 
Holstein model describing a single spin-degenerate elec- 
tronic state coupled to the vibration mode and to two su- 
perconducting electrodes with phase difference <f>. We fo- 
cus on the adiabatic regime, where the oscillator dynam- 
ics is slow on characteristic timescales of the electronic 
motion, and typically many oscillator quanta are excited 
under strong electron- vibration coupling. As shown be- 
low, this situation is analogous to a heavy Brownian par- 
ticle in a fast non-Ohmic fcrmionic environment The 



oscillator distribution function and the Josephson cur- 
rent can then be calculated by using a semiclassical de- 
scription of the oscillator dynamics. Thereby, a nonper- 
turbative treatment of the electron- vibration coupling is 
possible and controlled calculations can be performed in 
so far unexplored parameter regimes. Similar ideas have 
been employed before in the description of noncquilib- 
rium normal-state transport for this model^*^ which are 
here generalized to the superconducting case. We address 
the equilibrium case (no bias voltage), where the phase 
difference 4> is the relevant control parameter coupling to 
the oscillator's motion. Note that typical superconduct- 
ing gap scales are of the order of A w 1 meV, and in 
many experimentally studied cases? 3 - 2 - - — the relevant vi- 
brational energy scale is significantly smaller than A and 
our theory is directly applicable. 

The structure of this paper is as follows. In Sec. [Hi 
we discuss the model and introduce our semiclassical 
approach. In Sec. IIII1 we then consider the oscillator 
dynamics within the Born-Oppenheimer approximation, 
and we derive the Fokker-Planck equation for the distri- 
bution function in energy space. The approach is em- 
ployed to obtain the results presented in Sec. HVl fol- 
lowed by a discussion and some conclusions in Sec. fVl 
We mostly use units where e = h = ks = 1. 

II. MODEL AND SEMICLASSICAL APPROACH 
A. Model 

We start by describing a minimal model of a molec- 
ular quantum dot sandwiched by two superconducting 
leads. Similar to the normal-state case^ this model can 
capture essential aspects of the relevant physics in such 
devices. Writing the full Hamiltonian H = Hq + Ht + 
-f^ieads, the term Hq describes the isolated "molecule", 
including the vibration mode and its coupling to the 
electronic state. Ht refers to the electronic tunneling 
Hamiltonian connecting the molecular level to the super- 
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conducting electrodes, and -f/i ea ds describes the s-wave 
BCS superconducting leads with phase difference <f>. In 
this work, we only discuss the equilibrium case where 
both leads have the same chemical potential. Concern- 
ing Hq, we assume that only one spin-degenerate molec- 
ular electronic state is relevant, with (bare) energy eo. 
The corresponding fermion operator is d a for spin pro- 
jection a =t, I- The molecular dot is supposed to also 
host a dominant vibration mode of frequency f2, and we 
retain only this quantum oscillator mode with dimension- 
less canonically conjugate operators x and p. Following 
standard arguments^ usually the most important cou- 
pling (A) to d a is contained in a dot Hamiltonian of the 
form 



{dU<r-l)+^(p 2 +X 2 )- (1) 



In this representation, the oscillator acts as a charge 
parity detector, since it is displaced by the dot charge 
n = d\d a only for even n = {0,2}. It is convenient 
to employ the Nambu formalism, where fermion opera- 
tors are combined in the Nambu spinors d = (d^,d^) T 

and tpj k = (V'jfc.fi Yj(-k) \) T f° r electrons in the left or 
right lead (j = L/R) with momentum k. The leads are 
then described by a pair of BCS Hamiltonians, 



leads 



Y] ^nk ( e k a z + Acr *) V'jfei 



j,k 



(2) 



with normal-state dispersion tk and BCS gap A (taken 
real and positive) ; for simplicity, we consider identical su- 
perconductors. The standard Pauli matrices <J x , y , z and 
(To = diag(l, 1) act in Nambu space. Tunneling of elec- 
trons between the dot and the leads corresponds to 



Ht = 



E 

j=L,R=+, 



A,- 



t i>] k a z e ±Ur "f^d + h.c, (3) 



where we assume that both dot-to-lead tunnel couplings 
to are equal and /c-independent; the generalization to 
asymmetric cases is straightforward. The superconduct- 
ing phase difference </> enters via the phase factor dressing 
the tunnel matrix element. Finally, we define the hy- 
bridization energy T = 7rz^o|£o| 2 ; where vq = 2^ fc 5(efc) 
is the normal density of states in the leads. 

We next employ the real-time path integral 
techniquei^ to derive an effective action for the 
oscillator alone, i.e., we integrate out all electronic 
degrees of freedom. Although we study an equilibrium 
problem, it is technically easier to obtain the semiclas- 
sical limit from the real-time Keldysh technique^ We 
thus introduce the standard forward [backward] branch 
of the Keldysh contour, with oscillator trajectories 
xi(t) [x2{t)\- These define the classical trajectory 
x{t) = (xi + .T2)/2 and the quantum part y(t) = x\ — x%. 
The path-integral expression for the time evolution 
operator of the system then takes the form 



Z = / VxVy e 



i(S +S c ) 



(4) 



where the action of the uncoupled oscillator is 

S = - [ dt y{QT x x + fix), (5) 



and S e is an influence functional which results from trac- 
ing out all fermionic variables 



S e = -£Trln(<3 1 - ^cr z y 



(0) 



= -iTrlnG- 1 



(A/2)' 



Tr {Ga z y) n . 



The trace operation "Tr" extends over Nambu, Keldysh, 
and time (or energy) space, while the symbols "TVjv" 
( u Tik") used below will refer to a trace over Nambu 
(Keldysh) space only. G denotes the Keldysh Green's 
function (GF) of the dot for given classical trajectory 
(y = 0); the check notation (") indicates the 2x2 Keldysh 
structure. In terms of the Nambu spinors d, G(tx,t%) = 
— i(Tc[d(ti)d l! (t2)}) , where 7c is the time-ordering op- 
erator along the Keldysh contour. It is convenient to 
express the GF G(t 1 ,t 2 ) = &(t;r) with t = (h + t 2 )/2 
and t — ti — ti in the Wigncr ("mixed") representation. 
We will also frequently employ the Fourier transformed 
expression, G(t; r) = J du e^" 7 " G(t;u>). 

According to Eq. ((T|), for a given classical trajectory 
{x(t}} of the oscillator, the dot Keldysh GF can be ob- 
tained from the Dyson equation 



G~ 



Gn 



\x(t)a z 



(7) 



which formally represents an infinite-dimensional ma- 
trix equation in time (or energy) space and in Nambu- 
Keldysh space. The (inverse) GF in the absence of the 
electron-vibration coupling is 



G 1 = (id t - e a z )f z 



(8) 



where the Pauli matrices f x ,y,z a ct in Keldysh space. The 
self-energy S originates from the integration over the lead 
fermions. In frequency representation, the retarded and 
advanced components (in Nambu space) are given by^^ 

y/(u±i0+) 2 - A 2 

while the Keldysh component follows from the standard 
equilibrium relation, 

Z k (lj) = f(u) [Z R {u) - £ A (u;)] , f(u) = tanh(w/2T). 

(10) 

These components determine the Keldysh matrix struc- 
ture according to (y = ± denotes the upper/lower branch 
of the Keldysh contour) 

t vv ,{r) = \f |V~ T [vV R + v'V A + vv'Y*] ( W ). 

(H) 

Similarly, the Keldysh GF G can be decomposed into the 
retarded [G R ], advanced [G A ], and Keldysh [G K ] compo- 
nents. Note that so far our expressions are exact. 
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B. Semiclassical approach 

In this paper, our main interest concerns the adiabatic 
case of a slow oscillator, where f2 is the smallest physical 
frequency scale. In this limit, the kinetic term in Sq fa- 
vors small quantum fluctuations y(t), and a semiclassical 
approach expanding in {y(t)} becomes possible. (This 
approximation can also be justified in the limit of high 
temperatures.) For the normal case (A = 0), such an ap- 
proach has been worked out in detail for nonequilibrium 
transport in Refs. and |T(1 It constitutes a controlled 
appoximation for H < T and arbitrary A. In the super- 
conducting case, we instead require <C min(A, T). 

Within a semiclassical approach, we thus evaluate the 
action S e up to second order in the quantum amplitude 
y(t) while keeping the full nonlinear dependence on the 
classical trajectory x(t), 

S e = -iTi-lnG- 1 + + S {2) + 0(y 3 ). (12) 

The first-order term is S'e 1 '' = / dtT(t)y(t), with the to- 
tal force exerted by electrons on the oscillator 

T{t) = yTr w (G K (t,t)a z ) = A(l - (n(t))). (13) 
Using the Dyson equation (JSJ, some algebra yields 



T(t)=F e (t)- / dt'r)(t,t')x(t') 



(14) 



with the time-local part of the force, 

F e (t) = y Tr N [G$(t,t)a z ] + V (t,t)x(t). (15) 

Here the equal-time value of the damping kernel is 
i\ 2 f 

r){t,t) = — / dt'Tr NiK (G (t - t')f z a z G(t' ,t)a z ) . 



(16) 

The second term in Eq. (fbfj) describes retarded damp- 
ing, where the Wigner representation of the real-valued 
damping kernel rj(t,t') with t > t' is obtained by solving 
the equation 



1 \ , , A 2 v-^ fdu' 



(17) 



xTr 



N 



s=± 
K , 



A (t;u)' + su/2)cj z G K {t;uj' - suj/2)a 



-G$(t; J + soj/2)a z A(t; J - sw/2)a 2 



Here A = i(G R — G A ) is the full spectral function of the 
dot (including the electron- vibration coupling), while A$ 
refers to the corresponding A = case. The second-order 
noise term is from Eq. ((6|) given by 



S^ = - I dtdf y{t)K{t,t')y{t'), 



where the fluctuation kernel has the Wigner representa- 
tion 



= y / ^Tt N ([A + iG K ](t;cj' +uj/2) 



a z [A - iG K ] (t; (J - uj/2)a z S j . 



(19) 



Weak coupling limit: Fluctuation-dissipation 
theorem 



So far no approximations have been made in treating 
Si X ' 2 \ and the above expressions are exact. Before we 
address the adiabatic regime of small Q, it is instructive 
to briefly consider the case of small A but arbitrary il. In 
that case, the full GF G entering Eqs. (Tig]). ([T7|l and ([T9j) 
can be replaced by the free GF Go- Taking into account 
that G$(u) = /(w) [Gf(w)-G^(w)], cf. Eq. (HJ), we 
find from Eq. (fT7|| 



\2 p J. I 

x Trjy [A {lu' + u/2)a z A Q (uj' - lo/2)cf z ] . 



Some algebra shows that the fluctuation kernel K(uS) in 
Eq. (jTSJ) can also be expressed in terms of rj{uS), 



K{u) = uj \n B {uj) + 1] f?(w), 



(21) 



j/T 



1) 1 is the Bosc-Einstcin func- 



whcrc tib(w) = 
tion. 

Equation (T2T|) constitutes the well-known fluctuation- 
dissipation theorem^ for weak electron-vibration cou- 
pling and provides a consistency check for our formalism. 
In the normal state (A = 0), the damping kernel r)(uj) 
is often assumed to be a smooth function of w, which 
is then approximated by a constant, j]q = rj(uj = 0), 
according to Eq. (|20[) . Under this approximation, the 
damping constant entering the equation of motion is 
just ?7o/2, as follows from the resulting first-order ac- 
tion, = J dty[F e — (r]o/2)x]. In the high-temperature 
limit, the fluctuation kernel then describes white noise, 
K(t) = i]oT6(t). In the superconducting case (A ^ 0), 
however, the above kernels may not exhibit the assumed 
spectral smoothness. The presence of Andrcev bound 
states is known to cause singular behavior of the dot 
spectral function A(t;uj) in the subgap region |u>| < A. 
We will therefore take into account the electron damping 
and fluctuation effects on the oscillator's motion through- 
out the whole spectral range. Furthermore, we now go 
beyond the weak-coupling limit and consider a nonpcr- 
turbative theory in the electron-vibration coupling A. 
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III. ADIABATIC REGIME 
A. Born-Oppenheimer approximation 

Next we turn to the oscillator dynamics in the adi- 
abatic regime realized for Q -C min(r, E a ((j))), where 
E a (4>) is the Andreev level energy (see below). In partic- 
ular, we require Q <C A, which also implies that normal- 
state results do not follow from the expressions below by 
sending A — > 0. Since the oscillator dynamics is now 
much slower than the electronic motion, we can invoke 
the Born-Oppcnhcimcr (BO) approximation.— The dot 
GF G is thereby approximated by the adiabatic Green's 
function Q(t;uj) whose inverse for given t,tu follows from 
the Dyson equation 



g- 1 (t;u J ) = Go 1 (u J )-\cT z f z x(t), 



(22) 



which now is a simple 4x4 matrix (in Keldysh-Nambu 
space) relation. [We often denote Q(x(t);ui) = Q(t;ui).] 
This GF describes a noninteracting dot whose time- 
dependent energy level e(t) = eo + Xx(t) is determined by 
the instantaneous displacement x = x(t) of the oscillator. 
Equation (f2"2"j) is formally obtained from the Dyson equa- 
tion ([7]) for G{ti,t2) in the mixed representation (with 
t h2 =t±T/2), 



G{t-r) = G (r)+A / dt'Go(h-t')f z a z 
-t)x(t) + ...) 
-d t + ..) G(t;t'-t 2 ). 



x (x(t) + (t' - t)x(t) + ...) 
h - 1' , 



x 1 



Noting that d t corresponds to xd x , all derivative terms 
are of order 0(x) and should therefore be neglected 
within the BO approximation. Using the analogy to 
an effectively noninteracting quantum dot level, the re- 
tarded component of the adiabatic GF follows in the 
form^a 



R u)(l + a u ) + e(x)a z + a u Acos{(j)/2)a x 

V(x; Co) 

where we introduce the auxiliary quantity 

iT 



(23) 
(24) 



yJ{uj + iQ+Y - A 2 

and the denominator is given by 

V{x;lu) = w 2 (l + a u ) 2 -e 2 (x)-a 2 A 2 cos 2 (</>/2). (25) 

The resulting adiabatic spectral function A(x; u>) = A a + 
A c receives contributions from Andreev levels (A a , for 
|cj| < A) and from quasiparticle continuum states (A c , 
for \uj\ > A). The Andreev level spectral function is given 
by 



A a {x;uj) 



2tt 



-6 (\u\ - E a (x)) 



d L0 T>{x]u) 

x [lo(1 + a u ) + e(x)a z + a LO Acos(cj)/2)a x 



(26) 



where the ^-dependent Andreev level energy E a (x) £ 
[0, A) is a non-negative root of the equation T>(x; E a ) = 
0. The other components of the adiabatic GF then follow 
from Q A = [Q R Y and the equilibrium relation 



Q K {t;u) = -if{uj)A{x{t);u). 



(27) 



Next we address the corresponding adiabatic expres- 
sions for the damping and fluctuation kernels. Since the 
damping kernel in Eq. (|14p is multiplied by i, it is suffi- 
cient to replace G — > Q in the calculation of the damping 
kernel r\. Using Eqs. (TT()]) and (f2"2")l , the time- local force 
(|T5l) reads 



(28) 



while the damping [Eq. (|T7|) ] and fluctuation [Eq. (TT 
kernels are 

T?(t;w) = rj(x(t),0;u), (29) 
K{t;oj) = w[n B (oj) + 1] rj(x(t),x(t);uj), 

with a generalized "dissipation function" 

A 2 f rh / 

fj(x,x';u;) = ^- / ^—[f{J + u/2)-f(J-uo/2)\ 



Auj J 2ir 

^2 Tv n (A(x; J + slu/2)o- z A{x'; J - sw/2)o- z ) 

s=± 

= fj(x,x';—uj) = fj(x' ,x]uj). (30) 



For small A, the x-dependence in Eq. (|29[) can be ne- 
glected, and we recover the fluctuation-dissipation theo- 
rem (f2"Tj) . The decomposition A = A a + A c implies that 
V = Va + Vc and K = K a + K c separate into contributions 
from the Andreev level states and from continuum states. 
Mixed terms involving transitions between Andreev level 
and continuum states turn out to be always strongly sup- 
pressed in the adiabatic regime due to the presence of an 
energy threshold A — E a in the fcrmionic spectrum. As 
a result, such terms can safely be neglected. 

At this point, we pause and summarize what we have 
achieved so far. The total effective action of the oscillator 
within the BO approximation is 



S = - dty 



fl- 1 x-F(x)+ / dt' 7]{t,t')x{t') 



dtdt' y(t)K(t, t')y(t') - i Tr In £~ 



(31) 



where F(x) = — fix + F e {x) is the total potential force. 
The Wigncr representation of the dissipation and fluctu- 
ation kernels is 



r)(t;r) 
K(t;r) 



du 

— cos(ojt) r){x{t), 0; w), 

7T 

duo 



(32) 



— cos(cjt) uj coth(w/2T) fi(x(t), x{t); w), 
o 27r 
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with the generalized dissipation function fj in Eq. (|30[) . 

Before proceeding with the solution of the above 
stochastic problem, we briefly address the analytically 
tractable limit r 2> A. While the resulting expressions 
are not used in the full numerical solution in Sec. lIVl they 
are useful to develop intuition and to determine whether 
the undcrdamped vs ovcrdamped regime is realized, see 
Sec. IIII Dl For f> A, the Andreev level energy is! 



20 



E a (x) = Ay/l -7» sin 2 (0/2) 



(33) 



with the x-dependent effective transmission probability 

1 



T(x) 



1 + e 2 (x)/r 2 



(34) 



The Andreev level contribution to the fluctuation kernel 
K(x(t);uj) is then given by 



K a (x;u>) 



— S n {x)S(u-2nE a (x)). (35) 

n=0,±l 



With the Fermi function np(uj) = {e^l T + 1) 1 , we use 
the auxiliary functions 

So(z) = n F (E a (x))n F (-E a (x)){l-T(x))T 3 (x) 
4A 4 



■sin 4 (0/2), 



~ ±1 (x) = [n F ( T E a (x))] 2 T 3 (x) 



A 4 sin 2 4 
2El{x) 



The zero-frequency peak in K a (x;u>) is determined by 
quantum fluctuations of the Andreev level current^ 
when the reflectivity is finite, T(x) < 1. The contin- 
uum contribution to the above kernels exhibits only very 
weak dependence on uj and can be evaluated by tak- 
ing the uj — > limit. We find K c (x;uj) ~ Tj] c (x;u>) 
with t] c {x;uj) ~ (2A/r) 2 e _A / T . As a function of uj, 
both kernels r)(x(t);uj) and K(x(t);u>) exhibit a constant 
background due to the continuum states, responsible for 
Ohmic dissipation^ Superimposed on this Ohmic part, 
we have the Andreev level 5- type contributions. They 
include a peak at zero frequency. We then turn back to 
the full problem characterized by arbitrary ratio T/A. 



B. Fokker-Planck equation in energy space 

Following standard arguments^ we can transform the 
effective action pip to an equivalent Langevin equation. 
Using a Hubbard-Stratonovich transformation, 



-yKy/2 _ 



-tK- x £/2+i(Ly 



we introduce the auxiliary field £(f). Functional integra- 
tion over y then yields 



which enforces the Langevin equation 



C[x, £] = rr\i - F(x) + / dt' r)(t, t')x{t') - £(t) = 



(37) 

with Gaussian noise £(t). The stochastic noise field has 
zero mean and variance (£(t)£(t')) = K(t,t'). In what 
follows, we consider the weak damping limit,— where the 
damping force (oc x) is small compared to the potential 
force F{x). In this undcrdamped regime, the oscillator 
energy E varies slowly on the timescale Vt~ l , We quan- 
tify this condition in terms of the system parameters in 
Sec. lrTlDl 

We proceed by writing a Langevin equation for the 
slow energy variable E(t), which is averaged in time over 
the (energy-dependent) oscillator period Tg. Technically, 
by multiplying the Langevin equation by x and defining 
£ = x £, we find 

jt ( m + u{x) ) = ~f dt ' ^ t'ww^ + &)> 

(38) 

where U(x) is an effective oscillator potential with 
F{x) = -d x U{x), and = K(t,t')x{t)x{t'). 

The change of the oscillator energy 

E(t) = x 2 /(2tt) + U{x) 

is thus determined by the work done by the damping 
force —T]x and by the fluctuations. Averaging Eq. ([55)1 
over Te yields a Langevin equation in energy space, 

E(t) = -r,(E) + Z E (t), (€B(t)^(f)} = K(E)S(t - t'), 

(39) 

where rj(E) determines the energy dissipation rate and 
K(E) describes multiplicative (state-dependent) noise, 



n(E) 

K(E) 



Te dt r l 

— I dt' n(t,t')x(t)x(t'), (40) 



o t e 
Te dt 
Te 



dt' K(t,t')x(t)x{t'). 



The corresponding Fokker-Planck equation^ governing 
the energy distribution function w(E,t) of the oscillator 
is 

d t w{E,t) = d E L(E)w(E,t) + ^d E [K(E)w(E,t)} 

(41) 

A similar Fokker-Planck equation has been derived 
previously^i 1 ^ for the normal-state case. However, the 
kernels rj(t,t') and K(t,t') were approximated by time- 
local [oc 5(t — t')} expressions in Eq. (|40l) . 

The stationary solution of Eq. (|4Tj) is given by the gen- 
eralized Boltzmann distribution, 



Z = jvxV^e-^ K ~^l 2 Aet{g- l )8(z[x^]), (36) w(E) = jV K \E) exp (- J 



dE' 



T eS (E>) 



(42) 
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where T e g(E) = K{E)/2rj{E) is an effective energy- 
dependent temperature and M a normalization constant. 

In order to compute the fluctuation-dissipation co- 
efficients (|4"0"]) , we introduce a velocity- velocity corre- 
lation function at energy E. For a periodic solution 
x = x{f) = x(t + TE) of the undamped oscillator problem 
at given energy E, we take the correlator 

Q E (t; T ) = x(t + T/2)x(t- T /2)\ E (43) 

OO 

n=—oo 

where Q E , n = Q* E ,n = Qe.-u and Cl E = 2n/T E . Using 
Eq. ([311), we find 

C Te dt 

V(E) = / t=- y"QB,n(*)^(t),0;nn B ), (44) 

= / ^r£gU() (45) 

x nfig coth(nn B /2T) r/(x(i), x(t); nVt E ). 

While the above equations are straightforward to solve 
numerically in the case of a single-well potential, it is 
also possible to encounter bistable behavior as reported 
for the normal-state caseJ£ We will discuss the transi- 
tion from a single- well to a double well potential U(x) in 
detail in Sec. IIVI For the case of a double-well poten- 
tial U(x) with barrier height Eb, there are two solutions 
Wifi{E) defined within each well region (E < Eb), and 
a third solution u>3(E) applicable for energies above the 
barrier (E > Eb). These solutions have to be matched 
by boundary conditions'^ In particular, continuity im- 
poses wi(Eb) + W2(Eb) = w$(Eb), and the transition 
probability to each well at the separatrix should be equal, 
wi(E b ) = w 2 {E b ). 

C. Current 

In the adiabatic approximation, the Joscphson current 
is given by^ 

/doj 
— f(cj)a u Tr N (a x (g R (x;u;)) osc ) , 

(46) 

which involves time-averaging over an oscillator period 
T E for given E, followed by an average over the oscillator 
energy using the stationary distribution ()42[) . 



dEw(E) 



dt 
o t e 



(47) 



x 5 \m + u{ ' x) ~ E ) 

Analytic continuation then yields for the Joscphson cur- 
rent 

7(0) = -2TA 2 sin(0) £ a*, (V-\x;iv n )) wc (48) 

v n >0 



with fcrmion Matsubara frequencies v n = (2n+l)7rT (in- 
teger n). Equation (f24| yields at v = oti U = T/\/ A 2 + v 2 , 
and a similar result is obtained for V from Eq. ((231) . 



D. Underdamped regime 

In practice, the physically most relevant parameter 
regime corresponds to underdamped motion of the oscil- 
lator. This can be shown by an estimate for rj(E) given 
next. In Sec. IIVI we also show the full numerical result for 
r](E) to self-consistcntly verify that one indeed stays in 
the weak-damping limit. Our analytical estimates were 
obtained for T/A > 1. 

For given energy E, the Andreev level contribution 
rj a {E) is non-zero only if the oscillator path x = x E {t) 
passes through x = 0. We find 



Va(E) 

n 2 



(49) 



with the dimensionless number 



I Tp \ /T(O)Asin(0/2) 
g a = n F (E )n F {-E )j^ I — 



r 



where Eq = E a (0) denotes the bare Andreev level en- 
ergy. This estimate is obtained for zero Andreev level 
width 7a = and by neglecting the (.E-dependent) renor- 
malization of the oscillator frequency. (In the numerical 
analysis below, we use 7 a = 0.0IA.) 

On the other hand, the continuum contribution to the 
damping kernel rj is estimated by 



Vc{E) 

n 2 



TE 

riA 



-A/T 



(50) 



Note that the two contributions scale differently with E. 
The underdamped regime is realized when r](E)/n 2 < 1. 
It is straightforward to observe from the above expres- 
sions that for A < r, this condition is always fulfilled. 
The underdamped regime may cover even significantly 
larger electron- vibration couplings A. 



IV. RESULTS AND DISCUSSION 

Let us now describe results obtained from this semi- 
classical approach. We here only consider parameter 
sets consistent with the assumption of underdamped adi- 
abatic motion of the oscillator, see Sec. IIIIDl In ad- 
dition, we shall assume good coupling between dot and 
electrodes, T/A > 1, consistent with the fact that we ne- 
glect Coulomb interaction effects on the dot. Note that 
the opposite case T/A < 1 was studied in Ref. [Tol 

The numerical calculation goes as follows. We 
first compute the effective potential U(x) according to 
Eq. (12"5)) . Having determined the effective potential U (x), 
the calculation proceeds by computing Q E ^ n as defined in 
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FIG. 1: (Color online) Josephson current in units of eA/h 
vs superconducting phase difference <j> for the interacting case 
(A = 0.5A: blue circles) and for A = (dashed curve). The 
system parameters (in units of A) are F = 8, eo = —0.1, O = 
0.05, with temperature T = 0.2. Inset: Interaction correction 
to the current, <5/ p h = /(A) — J(A = 0), vs phase difference <j) 
for the data in the main panel. 




FIG. 3: (Color online) Energy distribution function (|42|l for 
<j> = (black dashed) and <j) — 0.8ir (blue solid curve), using 
the same parameter set as in Fig. [TJ The dotted curve shows a 
Boltzmann distribution for temperature T. Insets: Damping 
kernel r/(E) vs energy E for <j> = [top left] and for <j) = 0.8n 
[top right]. The corresponding effective temperature T e ff = 
K(E)/2r](E) is shown as a function of E for these two values 
of (j> in the bottom right inset: black filled circles are for <j> — 0, 
and blue open circles are for (j> = 0.8-7T. 




x 



FIG. 2: (Color online) Effective potential U(x) vs dimen- 
sionless oscillator coordinate x for = (black dashed) and 
4> = 0.8-7T (blue solid curve). System parameters are as in 
Fig- HI Inset: Andreev level spectrum vs x for the two quoted 
values of 4>. 



Eq. P5|) . This involves a numerical solution of the clas- 
sical equations of motion in the potential U(x), which 
are always periodic (the oscillation period Te is thereby 
obtained numerically). Subsequently we compute the 
damping kernel rj(E) using Eq. and the fluctuation 
kernel K(E) from Eq. ([45]) . These kernels then result in 
the probability distribution w(E) according to Eq. (p4"2"l) . 
and finally the Josephson current-phase relation is ob- 
tained from Eq. ([25]). 



A. Single-well case 

Figures [TJ [2] and [3] show our numerical results for the 
following set of system parameters: T = 8A, e = — 0.1A, 
£1 = 0.05A, A = 0.5A, with temperature put to T = 
0.2A. For this parameter set, we are in the weak-coupling 
(undcrdamped) regime, where the above formalism can 
be safely applied. The effective potential U(x) then has 
a single minimum for all values of the phase difference 
4>, see Fig. [2] for <f> — and <fi = 0.8ir. This single-well 
behavior of the effective oscillator potential surface can 
be rationalized by noting that the electron force F e (x) is 
here mainly determined by the continuum contribution, 
which in turn is almost insensitive to the phase difference 
</>. Interestingly, as seen in Fig. [TJ the Josephson current 
is basically not modified, with only a very small neg- 
ative interaction correction even for a relatively strong 
electron- vibration coupling A. The weak sensitivity of 
the current to A comes from a strong localization of the 
oscillator near the bottom of the effective potential U(x) 
at x = 0, see Fig. [2] 

The corresponding distribution functions w(E) for </> = 
and <f> = 0-8tt are shown in Fig. [3] The observed singu- 
lar behavior for small energies E is mainly determined by 
the factor K~ 1 (E) in Eq. (|12")) . For instance, for the con- 
tinuum contribution, one obtains K C (E) ~ 2Tr) c (E) oc 
E, and hence we find the scaling w(E) oc 1/E. The ap- 
proximately linear law K(E) oc E as E — > stays also 
valid when including the Andreev level contribution. In- 
deed, from Eq. 05]), we find K{E) m ETfj(x,x), with the 
average over phase space (at given energy E) defined as 

W^x) = fd 7^1r ] - Note that §dxp{x) « toE/Sl, 



8 




0.4- 



0.3 - 



0.2 - 



0.1 - 



//- r 



to 



OS- 



-0.1 - 



0.2 



0.4 0.6 
(|)/7t 



O.N 



\ 

\ 

\ 

'. \ 
v. \ 
■>. \ 
W 

: :\ 



FIG. 4: (Color online) Effective potential t/(a:) vs a; for = 
(black dashed), (f> — 0.65-7T (red dotted), and (f> = 0.9757T 
(blue solid curve). System parameters (with A = 1) are Y = 
4.8, e = -0.15, Q = 0.02, A = 0.4, T = 0.25. Inset: Andreev 
level spectrum vs x for these three values of 4>. 



FIG. 5: (Color online) Josephson current (in units of eA/h) 
vs (f> for the interacting case (A = 0.4A : blue circles) and 
for A = (black dashed curve). The thin dotted curve shows 
the result when x is held fixed at the global minimum of 
U(x). System parameters are as in Fig. [4] Inset: Interaction 
correction 8I p h vs phase difference. 



while fj(x, x) is only weakly dependent on E. We thus 
conclude again that K{E) oc E. 

For small <j), the main contribution to the oscillator 
damping comes from the continuum states, while for in- 
termediate <p the Andreev level contribution starts to 
dominate. This is explicitly seen in the two upper insets 
in Fig. |3j where rj(E) is shown for <fi = and <f> = 0.87T, 
respectively. For <f> = 0, we find a linear .E-dependence, 
which turns into a square-root dependence for = O.871", 
in accordance with Eqs. ([50)) and (|49[) . respectively. Fur- 
thermore, the bottom right inset of Fig. [3] shows that the 
effective temperature T e g(E) = K(E)/2r](E) is greatly 
enhanced for <j) = 0.87T due to Andreev-level current fluc- 
tuations. These fluctuations lead to stronger localization 
of the oscillator at low E. For = 0, we find T e s(E) ~ T, 
as expected when the continuum contributions dominate. 
Nevertheless, even then the oscillator distribution func- 
tion w{E) strongly deviates from the classical Boltzmann 
distribution of a free oscillator (shown in Fig. [3] for com- 
parison). This difference can be traced to the prcfactor 
K~ X (E) in Eq. 02]). 



B. Crossover to the double-well potential 

Let us next analyze a second parameter set, where we 
will encounter a nontrivial double-well behavior for the 
effective oscillator potential U(x). The transition from 
single- to double- well behavior is here induced by a vari- 
ation of the phase difference <$>, and one may therefore 
affect the conformational state of the molecule in a dis- 
sipationless manner in such a setup. The parameter set 
is given by fi = 0.02A, T = 4.8A, e = -0.15A, with 
electron-vibration coupling strength A = 0.4A. More- 
over, the temperature has been set to T — 0.25A. As 
illustrated in Fig. 21 we indeed find a transition between 



a single- and a double-well potential induced by a vari- 
ation of <p. Similar transitions (with associated bista- 
bilities) were reported for a two-level system instead of 
the oscillator ^ and for the noncquilibrium normal-state 
local Holstein models 

Although the continuum contribution to the electron 
force and thus to the effective potential U(x) still plays 
an overall dominant role, it is almost insensitive to vari- 
ations of cf>. The 0-tunable transition to a double-well 
potential shown in Figgis therefore caused by Andreev 
level contributions to F e {x). We note that the shape of 
U(x) is also sensitive to temperature through thermal 
occupation factors of the Andreev levels. The dynamical 
frequency Q e for the oscillator motion in the effective po- 
tential U (x) can be strongly renormalized away from the 
bare oscillator frequency Q. For the parameters in Fig. [4J 
we typically find £1e ~ 0.557. The cc-dependence of the 
Andreev level spectrum E a (x) in the adiabatic limit, i.e., 
with instantaneous x = x(t), is shown for several phases 
cf) in the inset of Fig. [4] Note that this spectrum is rather 
different from the featureless Andreev level spectrum for 
the first parameter set, see inset of Fig. 

Figure [5] shows the current-phase relation for this pa- 
rameter set. The Josephson current again exhibits an 
overall suppression due to the electron-vibration cou- 
pling as reported previouslyi^ ~ 30 i 39 The suppression is 
now more pronounced than in Fig. [U but still remains 
moderate. Moreover, the current-phase relation exhibits 
small yet characteristic cusps in the crossover region be- 
tween the single- and double-well situation (<j) w 0.6-7T 
to 0.77r), where Andreev level noise oc I 2 {(f) is strongly 
enhanced . 37 ' 38 However, the effect of switching between 
the two potential wells does not have a dramatic influence 
on the current-phase relation because the magnitude of 
the current is basically the same in each well: the coor- 



0.1 0.2 0.3 0.4 

E/A 

FIG. 6: (Color online) Energy distribution function (|42p 
(blue solid curve) for <f> — with parameters in Fig. [4] Note 
that the potential U(x) then has a single minimum. For 
comparison, the black dotted (black dashed) curve shows a 
Boltzmann (Bose-Einstein) distribution function for the given 
temperature T = 0.25A. The three insets show the corre- 
sponding energy-dependence of the kernels K(E) [top left], 
the damping kernel n{E) [top right], and the effective tem- 
perature T eS (E) = K(E)/2n(E) [bottom right]. 
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FIG. 7: (Color online) Same as Fig. [6] but for (j> = 0.65tt 
(crossover regime between single- and double- well potential). 
In the insets, results for the deeper and the shallower well 
below the potential barrier are shown by filled black circles 
and open red squares, respectively. 

dinates x\ t 2 of two local minima are almost symmetric 
with respect to i = 0. Indeed, we find Xi « — a; 2 and 
hence e(x±) « \x± « — 6(0:2) for strong coupling A and 
small eo- 

Figures |5J [7] and [5] then show the resulting energy dis- 
tribution function w(E) for the three values of the phase 
difference <j> considered in Fig. 31 respectively. Despite of 
the enhanced effective temperature T c g(E) as compared 
to T, we find again that w(E) strongly deviates from 
a classical (Boltzmann) distribution function. Interest- 
ingly, for energies E below the separatrix region, w{E) is 
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FIG. 8: (Color online) Same as Fig. but for <j> = 0.975tt, 
where U(x) in Fig. [4] corresponds to a double- well potential. 




FIG. 9: (Color online) Effective phase diagram in the A — eo 
plane for = 0.02 (units are such that A = 1). The main 
panel shows Acl(eo) for T = 0.1 (solid curves) and T = 0.4 
(dashed curves), where single- well behavior is found for A < 
A c i and double- well behavior starts to set in for A > A c i. 
The black (lower) curves are for T = 2, and the blue (upper) 
curves are for F = 8. Inset: Boundaries A c i (black lower 
curve) and A C 2 (red upper curve) for F = 2 and T = 0.1. For 
A > A C 2, the double- well behavior occurs for all values of the 
phase difference <f>. 



well approximated by an effective Bosc-Einstcin function. 
Indeed, we find that for J?«T, both w(E) and the Bose- 
Einstein function ub{E) obey the same equation. As a 
result, within this region, we find w(E) oc l/E instead of 
the Boltzmann dependence cx e~ E l T , implying a much 
stronger localization of the oscillator's motion near the 
bottom of the deeper well. Moreover, the damping n(E) 
is determined by both the continuum and Andreev level 
contributions, while the diffusion coefficient K(E) is es- 
sentially determined by the Andreev level contribution. 

Finally, we address the parameter regime where the 
described switching from single- to double-well behavior 
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in U (x) is found, see Fig. [9l For simplicity, we consider a 
fixed vibrational frequency, fl = 0.02 A. For given system 
parameters, when increasing A, we find from our numer- 
ical scheme that double-well behavior starts to appear 
at A = A c i for = jr. When further increasing A, the 
double-well behavior extends to a region with cf> < tt as 
well. A second scale A C 2 > A c i is then defined such that 
for A > A C 2, the double- well behavior is found for all <fi. In 
order to determine Ad^, it is therefore sufficient to probe 
for the single-to-double-well transition at the phase dif- 
ferences (f> = 7T and 4> = 0. The transition region between 
A c i and A C 2 is in fact rather narrow, cf . the inset of Fig. [§] 
Note that the location of the switching transition is quite 
sensitive to temperature. In particular, with increasing 
T, the boundary A c i,2 shifts to bigger A values. 

V. CONCLUSIONS 

We have shown that the adiabatic limit allows to make 
analytical progress for an important model of molecular 
electronics, the superconducting local Holstein model. 
It describes a spinless resonant electronic level coupled 
to a single boson (vibration mode), where the resonant 
level is coupled to two superconducting reservoirs with a 
phase difference <f>. The adiabatic limit is realized when 
the oscillator frequency f2 is smaller than both the su- 
perconducting gap A and the dot-to-lead hybridization 
energy scale T. This regime allows for a semiclassical 
Born-Oppcnhcimer-type treatment, where the electronic 
degrees of freedom can be integrated out and give rise 
to an effective oscillator potential U(x). Moreover, they 
cause dissipativc damping and a stochastic noise force. 
The most relevant parameter regime turns out to be the 
undcrdamped one, where it is appropriate to consider dif- 
fusion in energy space, and the effects of damping 
and noise [_ftT(i?)] can be taken into account within a stan- 
dard Fokkcr-Planck scheme. The resulting distribution 
function w(E) solving the Fokker-Planck equation can 
be obtained numerically with moderate effort, and al- 
lows us to obtain quantitative results within a controlled 
approximation for fi <C min(A,r). 

The method has been applied to a calculation of the 
Josephson current-phase relation I ((f) . While the result- 



ing corrections to the Josephson current are generally 
small in magnitude even for strong electron-vibration 
coupling A, they still cause some features when the effec- 
tive potential U (x) changes character. In particular, it is 
possible to induce a change from a single- to a double- well 
potential surface by variation of the phase difference (f>. 
Near the transition point, we predict enhanced Andrccv 
level noise. Similar transitions from single- to double- well 
effective potentials were also reported for the normal- 
state case^ but for a nonequilibrium situation where a 
finite bias voltage is applied and dissipation is unavoid- 
able. In our equilibrium case, the transition is induced 
by a variation of the superconducting phase and there- 
fore is dissipationlcss. However, the effect of switching 
from one to two minima in U[x) on the Josephson cur- 
rent I(4>) is much weaker here, which can be rationalized 
by noting that the two minima are located symmetri- 
cally and the Josephson current is essentially identical 
when the oscillator is close to a given minimum. It is 
also worth mentioning that the Josephson current seems 
always to be suppressed by the coupling to the oscillator, 
independent of the normal-state transmission probabil- 
ity through the junction, i.e., the interaction correction 
is negative, in contrast to what happens in the normal 
stated This conclusion has also been reached via per- 
turbation theory in the electron-vibration coupling A for 
the Josephson current ] 29 ' 39 

A quantity that is much more sensitive to the existence 
of two minima in the effective oscillator potential U(x) is 
the phonon distribution function. As we have discussed 
in Sec. HVl in the double-well case, the phonon distribu- 
tion function has a characteristic two-peak structure and 
displays strong phonon localization. It may be possible 
to access this quantity experimentally via resonant co- 
herent phonon spectroscopy techniques ) 40 ' 41 and thereby 
provide clear signatures of the predicted crossover from 
single- to double-well behavior. 
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